/***************************************************************************
 *   Copyright (C) 2008 by Chris Bornholdt,,,   *
 *   amelek@dell   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
 
#ifndef _CONVOLUTION_KERNELS_HPP_
#define _CONVOLUTION_KERNELS_HPP_

#include <cmath>

template<typename Element = double>
class rectangular_pulse
{
	public:
		typedef Element result_type;
		typedef double argument_type;
		
		rectangular_pulse(double min_domain = -1, double max_domain = 1, 
			double height = 1) : min_domain_(min_domain), max_domain_(max_domain),
			height_(height) {}
			
		double min_domain() const	{	return min_domain_;	}
		double max_domain() const	{	return max_domain_;	}
			
		result_type operator()(double x) const
		{
			if(min_domain_ <= x && x <= max_domain_)
				return height_;
			return 0;
		}

	private:
		double min_domain_, max_domain_, height_;
};

template<typename Element = double>
class triangular_pulse
{
	public:
		typedef Element result_type;
		typedef double argument_type;
		
		triangular_pulse(double center = 0, double width = 1, double height = 1) 
			: center_(center), width_(std::fabs(width)), height_(height) {}
			
		double min_domain() const	{	return center_ - width_;	}
		double max_domain() const	{	return center_ + width_;	}
			
		result_type operator()(double x) const
		{
			double delta = std::fabs(x - center_);
			if(delta > width_)
				return 0;
			return height_ * (width_ - delta) / width_;
		}

	private:
		double center_, width_, height_;
};

template<typename Element = double>
class gaussian_kernel
{
	static const double range_ = 3;
	public:
		typedef Element result_type;
		typedef double argument_type;
		
		gaussian_kernel(double sigma, double height) :
			sigma_(sigma), height_(height) {}
			
		double min_domain() const	{	return -range_;	}
		double max_domain() const	{	return range_;	}
		
		result_type operator()(double x) const
		{
			double tmp = -x * x / (2 * sigma_ * sigma_);
			return height_ * std::exp(tmp);
		}
		
	private:
		double sigma_, height_;	
};

template<typename Element = double>
class bicubic_kernel
{
	public:
		typedef Element result_type;
		typedef double argument_type;
		
		//or a = -0.75
		bicubic_kernel(double a = -0.50) : a_(a) {}
		
		double min_domain() const	{	return -2;	}
		double max_domain() const	{	return 2;	}
		
		result_type operator()(double x) const
		{
			x = std::fabs(x);
			//std::cout << x << std::endl;
			if(x <= 1)
				return (a_ + 2) * (x * x * x) - (a_ + 3) * (x * x) + 1;
			if(x <= 2)
				return a_ * x * x * x - 5 * a_ * x * x + 8 * a_ * x - 4 * a_;
			return 0;
		}
	
	private:
		double a_;
};

template<typename Element = double>
class lanczos_kernel
{
	public:
		typedef Element result_type;
		typedef double argument_type;
		
		lanczos_kernel(double a = 3) : a_(a) {}
		
		double min_domain() const	{	return -a_;	}
		double max_domain() const	{	return a_;	}
		
		result_type operator()(double x) const
		{
			static const double pi = 3.141528;
			if(x == 0)
				return 1;
			if(-a_ < x && x < a_)
				return a_ * sin(pi * x) * sin(pi * x / a_) / (pi * pi * x * x);
			return 0;
		}
		
	private:
		double a_;
};
				
#endif
